module ghost_interp_mod

  use spark_cubed_sphere
  use recon
  use params_mod
  use state_mod
  use diag_mod
  use static_mod
  use math_mod

  implicit none

  private

  public ghost_interp_init
  public ghost_interp
  public ghost_interp_topo
  public ghost_interp_topo_deriv
  public ghost_interp_final
  
  interface
    subroutine ghost_interp_interface(state, diag)
      import state_type
      import diag_type
      type(state_type), intent(inout) :: state
      type(diag_type ), intent(inout) :: diag
    end subroutine ghost_interp_interface
  end interface
  
  procedure(ghost_interp_interface), pointer :: ghost_interp => null()

  type poly_array_type
    type(recon_type), allocatable :: poly(:,:)
  end type poly_array_type

  type(poly_array_type), target :: gst_poly_arrays(6)
  real(r_kind), allocatable :: fo(:,:)

contains

  subroutine ghost_interp_init(mesh)

    type(cubed_sphere_mesh_type), intent(in) :: mesh

    integer idom, igs, ige, jgs, jge, pc, max_ngp
    integer i, j, p, is, ie, js, je, ierr
    real(r16) x0, y0, xc(sw), yc(sw), scalex, scaley

    call log_notice('Initialize polynomial interpolation for ghost points.', pid=proc%id)
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%initialized) then
        call mesh%get_params(idom=idom, igs=igs, ige=ige, jgs=jgs, jge=jge, pc=pc)
        allocate(gst_poly_arrays(idom)%poly(igs:ige,jgs:jge))
      
        associate (domain => mesh%domain(idom), dx => mesh%domain(idom)%dx, dy => mesh%domain(idom)%dy)
        scalex = 1.0_r16 / dx; scaley = 1.0_r16 / dy
        !$OMP PARALLEL DO PRIVATE(is,ie,js,je,x0,y0,xc,yc,p,ierr) COLLAPSE(2)
        do j = jgs, jge
          do i = igs, ige
            ! init procedure pointer
            gst_poly_arrays(idom)%poly(i,j)%init              => poly_square_init
            gst_poly_arrays(idom)%poly(i,j)%calc_recon_matrix => poly_square_calc_recon_matrix     
            gst_poly_arrays(idom)%poly(i,j)%reconstruct       => poly_square_reconstruct
            gst_poly_arrays(idom)%poly(i,j)%reconstruct_deriv => poly_square_reconstruct_deriv
            if (domain%ghost(i,j)%ngp > 0) then      
              is = i - rw; ie = i + rw; js = j - rw; je = j + rw
              if (is <  1) then; is =  1; ie =  1 + 2 * rw; end if
              if (js <  1) then; js =  1; je =  1 + 2 * rw; end if
              if (ie > nx) then; is = nx - 2 * rw; ie = nx; end if
              if (je > ny) then; js = ny - 2 * rw; je = ny; end if
              x0 = domain%x(pc,i,j); xc = domain%x(pc,is:ie,js)
              y0 = domain%y(pc,i,j); yc = domain%y(pc,is,js:je)
              ! Scale coordinates to [-rw,rw].
              xc = (xc - x0) * scalex; yc = (yc - y0) * scaley
              call gst_poly_arrays(idom)%poly(i,j)%init(nd=2, sw=sw, xc=xc, yc=yc, dx=real(dx,16), dy=real(dy,16), is=is, ie=ie, js=js, je=je)
              do p = 1, domain%ghost(i,j)%ngp
                call gst_poly_arrays(idom)%poly(i,j)%add_point(x=(domain%ghost(i,j)%x(p)-x0)*scalex, &
                                                               y=(domain%ghost(i,j)%y(p)-y0)*scaley)
              end do
              call gst_poly_arrays(idom)%poly(i,j)%calc_recon_matrix(ierr)
              if (ierr /= 0) call log_error('Failed to calculate ghost reconstruction matrix!')
            end if
          end do
        end do
        !$OMP END PARALLEL DO
        end associate
      end if
    end do

    call mesh%get_params(max_ngp=max_ngp)
    if (.not. allocated(fo)) allocate(fo(max_ngp,nVar))
    
    if(model_type==0)then
      ghost_interp => ghost_interp_adv
    elseif(model_type==1)then
      ghost_interp => ghost_interp_barotropic
    elseif(model_type==2)then
      ghost_interp => ghost_interp_baroclinic
    endif

  end subroutine ghost_interp_init
        
  subroutine ghost_interp_adv(state, diag)
    type(state_type), intent(inout) :: state
    type(diag_type ), intent(inout) :: diag

    integer idom, ingb, idom_ngb, ivar, ngp, ierr
    integer igs, ige, jgs, jge, irs, ire, jrs, jre
    integer ims, ime, jms, jme, ids, ide, jds, jde, kds, kde
    integer p, i, j, k, ps, is, js
    real(r_kind) uc, vc, us, vs, tmp(nVar)
    real(r_kind), pointer :: qQ(:,:,:,:,:)
    type(recon_type), pointer :: poly
    type(cubed_sphere_domain_type), pointer :: domain_ngb

    associate (mesh => state%mesh)
    ! Zero quadrature of halo cells.
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%initialized) then
        call mesh%get_params(idom=idom, ims=ims, ime=ime, jms=jms, jme=jme, ids=ids, ide=ide, jds=jds, jde=jde)
        associate (qQ => diag%domain(idom)%qQ)
        if (mesh%domain(idom)%touch_edges(1)) qQ(:,  ims:0  ,  jms:jme,:,:) = 0
        if (mesh%domain(idom)%touch_edges(2)) qQ(:,1+ide:ime,  jms:jme,:,:) = 0
        if (mesh%domain(idom)%touch_edges(3)) qQ(:,  ims:ime,  jms:0  ,:,:) = 0
        if (mesh%domain(idom)%touch_edges(4)) qQ(:,  ims:ime,1+jde:jme,:,:) = 0
        end associate
      end if
    end do
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%initialized) then
        call mesh%get_params(idom=idom, igs=igs, ige=ige, jgs=jgs, jge=jge, kds=kds, kde=kde)
        associate (domain => mesh%domain(idom),       &
                   ghost  => mesh%domain(idom)%ghost, &
                   q      => state%domain(idom)%q     )
        !$OMP PARALLEL DO PRIVATE(ngp, poly, ivar, ierr, fo, p, idom_ngb, ps, is, js, domain_ngb, qQ, tmp) COLLAPSE(3)
        do k = kds, kde
          do j = jgs, jge
            do i = igs, ige
              ngp = ghost(i,j)%ngp
              poly => gst_poly_arrays(idom)%poly(i,j)
              if (ngp == 0) cycle ! Skip cell without ghost points.
              do ivar = 1, nVar
#ifdef DEBUG
                call poly%reconstruct(q(poly%is:poly%ie,poly%js:poly%je,k:k,ivar), fo(1:ngp,ivar), ierr=ierr)
                if (ierr /= 0) call log_error('Failed to reconstruct ghost point!')
#else
                call poly%reconstruct(q(poly%is:poly%ie,poly%js:poly%je,k:k,ivar), fo(1:ngp,ivar))
#endif                
              end do
              ! Copy ghost points to quadrature points.
              do p = 1, ngp
                idom_ngb = ghost(i,j)%idom(p)
                ps       = ghost(i,j)%ps  (p)
                is       = ghost(i,j)%is  (p)
                js       = ghost(i,j)%js  (p)
                ! Project velocity vector from this domain coordinate system to neightbor domain.
                domain_ngb => mesh%domain(idom_ngb); qQ => diag%domain(idom_ngb)%qQ
                tmp(1) = fo(p,1) / mesh%get_gst_Jh(idom,p,i,j,k) * mesh%get_dom_Jh(idom_ngb,ps,is,js,k)
                qQ(ps,is,js,k,:) = qQ(ps,is,js,k,:) + ghost(i,j)%wgt(p) * tmp
              end do
            end do
          end do
        end do
        !$OMP END PARALLEL DO
        end associate
      end if
    end do
    ! Calculate cell integrated values.
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%initialized) then
        call mesh%get_params(idom=idom, irs=irs, ire=ire, jrs=jrs, jre=jre, kds=kds, kde=kde)
        associate (domain => mesh%domain(idom),   &
                   q      => state%domain(idom)%q,& 
                   qQ     => diag%domain(idom)%qQ)
        !$OMP PARALLEL DO COLLAPSE(4)
        do ivar = 1, nVar
          do k = kds, kde
            do j = jrs, jre
              do i = irs, ire
                if (domain%cell_type(i,j) /= domain_cell) then
                  q(i,j,k,ivar) = domain%cell_quad(qQ(:,i,j,k,ivar))
                end if
              end do
            end do
          end do
        end do
        !$OMP END PARALLEL DO
        end associate
      end if
    end do
    end associate

  end subroutine ghost_interp_adv

  subroutine ghost_interp_barotropic(state, diag)
    type(state_type), intent(inout) :: state
    type(diag_type ), intent(inout) :: diag

    integer idom, ingb, idom_ngb, ivar, ngp, ierr
    integer igs, ige, jgs, jge, irs, ire, jrs, jre
    integer ims, ime, jms, jme, ids, ide, jds, jde, kds, kde
    integer p, i, j, k, ps, is, js
    real(r_kind) uc, vc, us, vs, tmp(nVar)
    real(r_kind), pointer :: qQ(:,:,:,:,:)
    type(recon_type), pointer :: poly
    type(cubed_sphere_domain_type), pointer :: domain_ngb

    associate (mesh => state%mesh)
    ! Zero quadrature of halo cells.
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%initialized) then
        call mesh%get_params(idom=idom, ims=ims, ime=ime, jms=jms, jme=jme, ids=ids, ide=ide, jds=jds, jde=jde)
        associate (qQ => diag%domain(idom)%qQ)
        if (mesh%domain(idom)%touch_edges(1)) qQ(:,  ims:0  ,  jms:jme,:,:) = 0
        if (mesh%domain(idom)%touch_edges(2)) qQ(:,1+ide:ime,  jms:jme,:,:) = 0
        if (mesh%domain(idom)%touch_edges(3)) qQ(:,  ims:ime,  jms:0  ,:,:) = 0
        if (mesh%domain(idom)%touch_edges(4)) qQ(:,  ims:ime,1+jde:jme,:,:) = 0
        end associate
      end if
    end do
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%initialized) then
        call mesh%get_params(idom=idom, igs=igs, ige=ige, jgs=jgs, jge=jge, kds=kds, kde=kde)
        associate (domain => mesh%domain(idom),       &
                   ghost  => mesh%domain(idom)%ghost, &
                   q      => state%domain(idom)%q     )
        !$OMP PARALLEL DO PRIVATE(ngp, poly, ivar, ierr, fo, p, idom_ngb, ps, is, js, uc, vc, us, vs, domain_ngb, qQ, tmp) COLLAPSE(3)
        do k = kds, kde
          do j = jgs, jge
            do i = igs, ige
              ngp = ghost(i,j)%ngp
              poly => gst_poly_arrays(idom)%poly(i,j)
              if (ngp == 0) cycle ! Skip cell without ghost points.
              do ivar = 1, nVar
#ifdef DEBUG
                call poly%reconstruct(q(poly%is:poly%ie,poly%js:poly%je,k:k,ivar), fo(1:ngp,ivar), ierr=ierr)
                if (ierr /= 0) call log_error('Failed to reconstruct ghost point!')
#else
                call poly%reconstruct(q(poly%is:poly%ie,poly%js:poly%je,k:k,ivar), fo(1:ngp,ivar))
#endif                
              end do
              ! Copy ghost points to quadrature points.
              do p = 1, ngp
                idom_ngb = ghost(i,j)%idom(p)
                ps       = ghost(i,j)%ps  (p)
                is       = ghost(i,j)%is  (p)
                js       = ghost(i,j)%js  (p)
                ! Project velocity vector from this domain coordinate system to neightbor domain.
                domain_ngb => mesh%domain(idom_ngb); qQ => diag%domain(idom_ngb)%qQ
                uc = fo(p,2) / fo(p,1); vc = fo(p,3) / fo(p,1)
                call contrav_to_spherev(mesh%get_gst_Ah (idom    ,p ,i ,j ,k), uc, vc, us, vs)
                call spherev_to_contrav(mesh%get_dom_iAh(idom_ngb,ps,is,js,k), us, vs, uc, vc)
                tmp(1) = fo(p,1) / mesh%get_gst_Jh(idom,p,i,j,k) * mesh%get_dom_Jh(idom_ngb,ps,is,js,k)
                tmp(2) = tmp(1) * uc
                tmp(3) = tmp(1) * vc
                qQ(ps,is,js,k,:) = qQ(ps,is,js,k,:) + ghost(i,j)%wgt(p) * tmp
              end do
            end do
          end do
        end do
        !$OMP END PARALLEL DO
        end associate
      end if
    end do
    ! Calculate cell integrated values.
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%initialized) then
        call mesh%get_params(idom=idom, irs=irs, ire=ire, jrs=jrs, jre=jre, kds=kds, kde=kde)
        associate (domain => mesh%domain(idom),   &
                   q      => state%domain(idom)%q,& 
                   qQ     => diag%domain(idom)%qQ)
        !$OMP PARALLEL DO COLLAPSE(4)
        do ivar = 1, nVar
          do k = kds, kde
            do j = jrs, jre
              do i = irs, ire
                if (domain%cell_type(i,j) /= domain_cell) then
                  q(i,j,k,ivar) = domain%cell_quad(qQ(:,i,j,k,ivar))
                end if
              end do
            end do
          end do
        end do
        !$OMP END PARALLEL DO
        end associate
      end if
    end do
    end associate

  end subroutine ghost_interp_barotropic

  subroutine ghost_interp_baroclinic(state, diag)
    type(state_type), intent(inout) :: state
    type(diag_type ), intent(inout) :: diag

    integer idom, ingb, idom_ngb, ivar, ngp, ierr
    integer igs, ige, jgs, jge, irs, ire, jrs, jre
    integer ims, ime, jms, jme, ids, ide, jds, jde, kds, kde
    integer pc, pqs, pqe
    integer p, i, j, k, ps, is, js
    real(r_kind) uc, vc, us, vs, w, rhod, pt, mr, sqrtGrho, tmp(nVar)
    real(r_kind), pointer :: qQ(:,:,:,:,:)
    type(recon_type), pointer :: poly
    type(cubed_sphere_domain_type), pointer :: domain_ngb

    associate (mesh => state%mesh)
    ! Zero quadrature of halo cells.
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%initialized) then
        call mesh%get_params(idom=idom, ims=ims, ime=ime, jms=jms, jme=jme, ids=ids, ide=ide, jds=jds, jde=jde)
        associate (qQ => diag%domain(idom)%qQ)
        if (mesh%domain(idom)%touch_edges(1)) qQ(:,  ims:0  ,  jms:jme,:,:) = 0
        if (mesh%domain(idom)%touch_edges(2)) qQ(:,1+ide:ime,  jms:jme,:,:) = 0
        if (mesh%domain(idom)%touch_edges(3)) qQ(:,  ims:ime,  jms:0  ,:,:) = 0
        if (mesh%domain(idom)%touch_edges(4)) qQ(:,  ims:ime,1+jde:jme,:,:) = 0
        end associate
      end if
    end do
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%initialized) then
        call mesh%get_params(idom=idom, igs=igs, ige=ige, jgs=jgs, jge=jge, kds=kds, kde=kde, pc=pc)
        associate (domain => mesh%domain(idom),       &
                   Jab    => mesh%domain(idom)%Jab,   &
                   ghost  => mesh%domain(idom)%ghost, &
                   q      => state%domain(idom)%q     )
        !$OMP PARALLEL DO PRIVATE(ngp, poly, ivar, ierr, fo, p, idom_ngb, ps, is, js, &
        !$OMP uc, vc, us, vs, w, rhod, pt, mr, sqrtGrho, domain_ngb, qQ, tmp) COLLAPSE(3)
        do k = kds, kde
          do j = jgs, jge
            do i = igs, ige
              ngp = ghost(i,j)%ngp
              poly => gst_poly_arrays(idom)%poly(i,j)
              if (ngp == 0) cycle ! Skip cell without ghost points.
              do ivar = 1, nVar
#ifdef DEBUG
                call poly%reconstruct(q(poly%is:poly%ie,poly%js:poly%je,k:k,ivar), fo(1:ngp,ivar), ierr=ierr)
                if (ierr /= 0) call log_error('Failed to reconstruct ghost point!')
#else
                call poly%reconstruct(q(poly%is:poly%ie,poly%js:poly%je,k:k,ivar), fo(1:ngp,ivar))
#endif                
              end do
              ! Copy ghost points to quadrature points.
              do p = 1, ngp
                idom_ngb = ghost(i,j)%idom(p)
                ps       = ghost(i,j)%ps  (p)
                is       = ghost(i,j)%is  (p)
                js       = ghost(i,j)%js  (p)
                ! Project velocity vector from this domain coordinate system to neightbor domain.
                domain_ngb => mesh%domain(idom_ngb)
                qQ => diag%domain(idom_ngb)%qQ
                sqrtGrho = fo(p,1) + fo(p,6)

#ifdef DEBUG
                if(sqrtGrho==0)then
                  print*,proc%id,idom,p,i,j,k
                  do iVar = 1, nVar
                    print*,q(poly%is:poly%ie,poly%js:poly%je,k:k,ivar)
                  enddo
                  do iVar = 1, nVar
                    print*,q(poly%is:poly%ie,poly%js:poly%je,k-1:k-1,ivar)
                  enddo
                  print*,''
                  print*,fo(p,:)
                  print*, 'Mass is 0 during ghost interpolation ps, is, js, i, j, k ', ps, is, js, i, j, k
                  stop
                endif
#endif
                
                rhod = fo(p,1)
                uc   = fo(p,2) / sqrtGrho
                vc   = fo(p,3) / sqrtGrho
                w    = fo(p,4) / sqrtGrho
                pt   = fo(p,5) / fo(p,1)
                mr   = fo(p,6) / fo(p,1)
                call contrav_to_spherev(mesh%get_gst_Ah (idom    ,p ,i ,j ,k), uc, vc, us, vs)
                call spherev_to_contrav(mesh%get_dom_iAh(idom_ngb,ps,is,js,k), us, vs, uc, vc)
                tmp(1)   = rhod
                tmp(5)   = tmp(1) * pt
                tmp(6)   = tmp(1) * mr
                sqrtGrho = tmp(1) + tmp(6)
                tmp(2)   = sqrtGrho * uc
                tmp(3)   = sqrtGrho * vc
                tmp(4)   = sqrtGrho * w
                qQ(ps,is,js,k,:) = qQ(ps,is,js,k,:) + ghost(i,j)%wgt(p) * tmp
              end do
            end do
          end do
        end do
        !$OMP END PARALLEL DO
        end associate
      end if
    end do
    ! Calculate cell integrated values.
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%initialized) then
        call mesh%get_params(idom=idom, irs=irs, ire=ire, jrs=jrs, jre=jre, kds=kds, kde=kde, pc=pc, pqs=pqs, pqe=pqe)
        associate (domain => mesh%domain(idom),    &
                   Jab    => mesh%domain(idom)%Jab,&
                   q      => state%domain(idom)%q, & 
                   qQ     => diag%domain(idom)%qQ  )
        !$OMP PARALLEL DO COLLAPSE(4)
        do ivar = 1, nVar
          do k = kds, kde
            do j = jrs, jre
              do i = irs, ire
                if (domain%cell_type(i,j) /= domain_cell) then
                  q(i,j,k,ivar) = domain%cell_quad( Jab(pqs:pqe,i,j,k) * qQ(:,i,j,k,ivar) ) / Jab(pc,i,j,k)
                end if
              end do
            end do
          end do
        end do
        !$OMP END PARALLEL DO
        end associate
      end if
    end do
    end associate

  end subroutine ghost_interp_baroclinic
  
  subroutine ghost_interp_topo(static)

    type(static_type), intent(inout) :: static
    
    integer idom, ingb, idom_ngb, ivar, ngp, max_ngp, ierr
    integer igs, ige, jgs, jge, irs, ire, jrs, jre
    integer ims, ime, jms, jme, ids, ide, jds, jde, kds, kde
    integer p, i, j, k, ps, is, js
    real(r_kind) uc, vc, us, vs, tmp(3)
    real(r_kind), pointer :: zsQ(:,:,:)
    real(r_kind), pointer :: zsG(:)
    
    type(recon_type), pointer :: poly
    
    real(r_kind), allocatable :: fi(:,:,:)
    real(r_kind), allocatable :: fo(:)
    
    !call log_notice('Calculating topo ghost interpolation.', pid=proc%id)

    associate (mesh => static%mesh)
      call mesh%get_params(max_ngp=max_ngp)
      if (.not. allocated(fi)) allocate(fi(-rw:rw,-rw:rw,1))
      if (.not. allocated(fo)) allocate(fo(max_ngp))
      
      ! Zero quadrature of halo cells.
      do idom = 1, mesh%ndom
        if (mesh%domain(idom)%initialized) then
          call mesh%get_params(idom=idom, ims=ims, ime=ime, jms=jms, jme=jme, ids=ids, ide=ide, jds=jds, jde=jde)
          associate (zsQ => static%domain(idom)%zsQ)
          if (mesh%domain(idom)%touch_edges(1)) zsQ(:,  ims:0  ,  jms:jme) = 0
          if (mesh%domain(idom)%touch_edges(2)) zsQ(:,1+ide:ime,  jms:jme) = 0
          if (mesh%domain(idom)%touch_edges(3)) zsQ(:,  ims:ime,  jms:0  ) = 0
          if (mesh%domain(idom)%touch_edges(4)) zsQ(:,  ims:ime,1+jde:jme) = 0
          end associate
        end if
      end do
      do idom = 1, mesh%ndom
        if (mesh%domain(idom)%initialized) then
          call mesh%get_params(idom=idom, igs=igs, ige=ige, jgs=jgs, jge=jge, kds=kds, kde=kde)
          associate (domain => mesh%domain(idom),       &
                     ghost  => mesh%domain(idom)%ghost, &
                     zs     => static%domain(idom)%zs)
          !$OMP PARALLEL DO PRIVATE(ngp, poly, fi, fo, p, idom_ngb, ps, is, js, zsQ, zsG) COLLAPSE(2)
          do j = jgs, jge
            do i = igs, ige
              ngp = ghost(i,j)%ngp
              if (ngp == 0) cycle ! Skip cell without ghost points.
              poly => gst_poly_arrays(idom)%poly(i,j)
              fi(:,:,1) = zs(poly%is:poly%ie,poly%js:poly%je)
              call poly%reconstruct(fi(:,:,1:1), fo(1:ngp))
              ! Copy ghost points to quadrature points.
              do p = 1, ngp
                idom_ngb = ghost(i,j)%idom(p)
                ps       = ghost(i,j)%ps  (p)
                is       = ghost(i,j)%is  (p)
                js       = ghost(i,j)%js  (p)
                ! Project velocity vector from this domain coordinate system to neightbor domain.
                zsQ  => static%domain(idom_ngb)%zsQ
                zsG  => ghost(i,j)%zs
                zsG(p) = fo(p)
                !zsG(p) = zs(i,j)
                zsQ(ps,is,js) = zsG(p)
              end do
            end do
          end do
          !$OMP END PARALLEL DO
          end associate
        end if
      end do
      
      ! Calculate cell integrated values.
      do idom = 1, mesh%ndom
        if (mesh%domain(idom)%initialized) then
          call mesh%get_params(idom=idom, irs=irs, ire=ire, jrs=jrs, jre=jre, kds=kds, kde=kde)
          associate (domain => mesh%domain(idom),      &
                     zs     => static%domain(idom)%zs, &
                     zsQ    => static%domain(idom)%zsQ )
          !$OMP PARALLEL DO COLLAPSE(4)
          do ivar = 1, nVar
            do k = kds, kde
              do j = jrs, jre
                do i = irs, ire
                  if (domain%cell_type(i,j) /= domain_cell) then
                    zs(i,j) = domain%cell_quad(zsQ(:,i,j))
                  end if
                end do
              end do
            end do
          end do
          !$OMP END PARALLEL DO
          end associate
        end if
      end do
    end associate
    
  end subroutine ghost_interp_topo
  
  subroutine ghost_interp_topo_deriv(static)

    type(static_type), intent(inout) :: static
    
    integer idom, ingb, idom_ngb, ivar, ngp, max_ngp, ierr
    integer igs, ige, jgs, jge, irs, ire, jrs, jre
    integer ims, ime, jms, jme, ids, ide, jds, jde, kds, kde
    integer p, i, j, k, ps, is, js
    real(r_kind) uc, vc, us, vs, tmp(3)
    real(r_kind), pointer :: dzsdxQ(:,:,:)
    real(r_kind), pointer :: dzsdxG(:)
    real(r_kind), pointer :: dzsdyQ(:,:,:)
    real(r_kind), pointer :: dzsdyG(:)
    type(recon_type), pointer :: poly
    
    real(r_kind) dzsdlon, dzsdlat
    real(r_kind) Ah(2,2), iAh(2,2)
    
    real(r_kind), allocatable :: fi(:,:,:)
    real(r_kind), allocatable :: fx(:)
    real(r_kind), allocatable :: fy(:)
    
    call log_notice('Calculating topo deriv ghost interpolation.', pid=proc%id)

    associate (mesh => static%mesh)
      call mesh%get_params(max_ngp=max_ngp)
      if (.not. allocated(fi)) allocate(fi(-rw:rw,-rw:rw,1))
      if (.not. allocated(fx)) allocate(fx(max_ngp))
      if (.not. allocated(fy)) allocate(fy(max_ngp))
      
      ! Zero quadrature of halo cells.
      do idom = 1, mesh%ndom
        if (mesh%domain(idom)%initialized) then
          call mesh%get_params(idom=idom, ims=ims, ime=ime, jms=jms, jme=jme, ids=ids, ide=ide, jds=jds, jde=jde)
          associate (dzsdxQ => static%domain(idom)%dzsdxQ)
          if (mesh%domain(idom)%touch_edges(1)) dzsdxQ(:,  ims:0  ,  jms:jme) = 0
          if (mesh%domain(idom)%touch_edges(2)) dzsdxQ(:,1+ide:ime,  jms:jme) = 0
          if (mesh%domain(idom)%touch_edges(3)) dzsdxQ(:,  ims:ime,  jms:0  ) = 0
          if (mesh%domain(idom)%touch_edges(4)) dzsdxQ(:,  ims:ime,1+jde:jme) = 0
          end associate
          associate (dzsdyQ => static%domain(idom)%dzsdyQ)
          if (mesh%domain(idom)%touch_edges(1)) dzsdyQ(:,  ims:0  ,  jms:jme) = 0
          if (mesh%domain(idom)%touch_edges(2)) dzsdyQ(:,1+ide:ime,  jms:jme) = 0
          if (mesh%domain(idom)%touch_edges(3)) dzsdyQ(:,  ims:ime,  jms:0  ) = 0
          if (mesh%domain(idom)%touch_edges(4)) dzsdyQ(:,  ims:ime,1+jde:jme) = 0
          end associate
        end if
      end do
      do idom = 1, mesh%ndom
        if (mesh%domain(idom)%initialized) then
          call mesh%get_params(idom=idom, igs=igs, ige=ige, jgs=jgs, jge=jge, kds=kds, kde=kde)
          associate (domain => mesh%domain(idom),       &
                     ghost  => mesh%domain(idom)%ghost, &
                     zs     => static%domain(idom)%zs   )
          k = 1
          !$OMP PARALLEL DO PRIVATE(ngp, poly, fi, fx, fy, p, idom_ngb, ps, is, js, &
          !$OMP                     dzsdlon, dzsdlat, dzsdxQ, dzsdxG, dzsdyQ, dzsdyG, iAh, Ah) COLLAPSE(2)
          do j = jgs, jge
            do i = igs, ige
              ngp = ghost(i,j)%ngp
              poly => gst_poly_arrays(idom)%poly(i,j)
              if (ngp == 0) cycle ! Skip cell without ghost points.
              fi(:,:,1) = zs(poly%is:poly%ie,poly%js:poly%je)
              call poly%reconstruct_deriv(fi(:,:,1:1), fx(1:ngp), fy(1:ngp))
              ! Copy ghost points to quadrature points.
              do p = 1, ngp
                idom_ngb = ghost(i,j)%idom(p)
                ps       = ghost(i,j)%ps  (p)
                is       = ghost(i,j)%is  (p)
                js       = ghost(i,j)%js  (p)
                ! Project velocity vector from this domain coordinate system to neightbor domain.
                dzsdxQ => static%domain(idom_ngb)%dzsdxQ
                dzsdyQ => static%domain(idom_ngb)%dzsdyQ
                dzsdxG => ghost(i,j)%dzsdx
                dzsdyG => ghost(i,j)%dzsdy

                dzsdxG(p) = fx(p)
                dzsdyG(p) = fy(p)
                
                ! Unify topo derivative
                iAh = mesh%get_gst_iAh(idom,p,i,j,k)
                call convert_hor_deriv_cube_to_sph(dzsdlon,dzsdlat,dzsdxG(p),dzsdyG(p),iAh)
                
                Ah = mesh%get_dom_Ah(idom_ngb,ps,is,js,k)
                call convert_hor_deriv_sph_to_cube(dzsdxQ(ps,is,js),dzsdyQ(ps,is,js),dzsdlon,dzsdlat,Ah)
              end do
            end do
          end do
          !$OMP END PARALLEL DO
          end associate
        end if
      end do
    end associate
  end subroutine ghost_interp_topo_deriv

  subroutine ghost_interp_final()

    integer idom

    do idom = 1, mesh%ndom
      if (allocated(gst_poly_arrays(idom)%poly)) deallocate(gst_poly_arrays(idom)%poly)
    end do
    if (allocated(fo)) deallocate(fo)

  end subroutine ghost_interp_final

end module ghost_interp_mod
